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Abstract 

In this paper, we propose an incremental algorithm for computing cylindrical al- 
gebraic decompositions. The algorithm consists of two parts: computing a complex 
cylindrical tree and refining this complex tree into a cylindrical tree in real space. The 
incrementality comes from the first part of the algorithm, where a complex cylindrical 
tree is constructed by refining a previous complex cylindrical tree with a polynomial 
constraint. We have implemented our algorithm in Maple. The experimentation shows 
that the proposed algorithm outperforms existing ones for many examples taken from 
the literature. 

1 Introduction 

Cylindrical algebraic decomposition (CAD) is a fundamental tool in real algebraic geom- 
etry. It was invented by G.E. Collins in 1973 jT5] for solving real quantifier elimination 
(QE) problems. In the last forty years, following Collins' original projection-lifting scheme, 
many enhancements have been performed in order to ameliorate the efficiency of CAD 
construction, including adjacency and clustering techniques [T], improved projection meth- 
ods [MlllSlllll], partially buih CADs [HlllTllS^, improved stack construction [18], efla- 
cient projection orders pO], making use of equational constraints [El [291 H [30] , and so on. 
Moreover, CADs can be computed by several software packages, such as Qepcad [^ [5]. 
Mathematica [Ml [31], Redlog [H] and SyNRAC [35]. 

In [TJ], together with B. Xia and L. Yang, we presented a different way for computing 
CADs based on triangular decomposition of polynomial systems. In that paper, we intro- 
duced the concept of cylindrical decomposition of the complex space (CCD), from which a 
CAD can be easily derived. The concept of CCD is reviewed in Section [2] In the rest of 
the present paper, we use TCAD to denote CAD based on triangular decompositions while 
PCAD refers to CAD based on Collins' projection-lifting scheme. 

The CCD part of TCAD can be seen as an enhanced projection phase of PCAD. 
However, w.r.t. PCAD (especially when the projection operator is using Collins' [T5| or 
Hong's [23]), the "case discussion" scheme of TCAD avoids unnecessary computations that 
projection operator performs on unrelated branches. In addition, one observes that the 
reason why McCallum's [28] (including Brown's [1]) projection operators may fail for some 
examples is due to the fact that they are missing a "case discussion" scheme. McCallum's 
operator relies on the assumption that generically all coefficients of a polynomiaj^ will not 
vanish simultaneously above a positive-dimensional component. If this assumption fails, 
then this operator is replaced by Collins-Hong projection-operator [53]. The fact that all 
coefficients of polynomial could vanish simultaneously above some component is never a 
problem in TCAD. For this reason, we view it as an improvement of previous works. 

^Morc precisely, a multivariate polynomial regarded as a univariate one with respect to its main variable. 
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Trying to use sophisticated algebraic elimination techniques to improve CAD construc- 
tions is not a new idea. In papers [51 [35] , the authors investigated how to use Grobner 
bases to preprocess the input system in order to make the subsequent CAD computations 
more efficient. The main difference between these two works and the work of |14| is that 
the former approach is about preprocessing input for CAD while the latter one presents a 
different way of constructing CADs. 

In |14j , the focus was on how to apply triangular decomposition techniques to compute 
CADs. To this end, lots of existing high-level routines were used to facilitate explaining 
ideas. These high-level routines involve many black-boxes, which hide many unnecessary 
or redundant computations. As a result, the computation time of TCAD is much higher 
than that of PCAD, although TCAD computes usually less cells [lO] . 

In the present paper, we abandon those black-boxes and compute TCAD from scratch. 
It turns out that the key solution for avoiding redundant computations is to compute CCD 
in an incremental manner. The same motivation and a similar strategy appeared in [321 112] 
in the context of triangular decomposition of algebraic sets. The core operation of such an 
incremental algorithm is an Intersect operation, which refines an existing cylindrical tree 
w.r.t. a polynomial. We dedicate Section |4] to presenting a complete incremental algorithm 
for computing TCAD by means of this Intersect operation. 

In |35| . the author presented an algorithm for computing with semi-algebraic sets rep- 
resented by cylindrical algebraic formulas. That algorithm also allows computing CAD in 
an incremental manner. The underlying technique is based on the projection-lifting scheme 
where one first computes projection factor sets by a global projection operator. In contrast, 
the incremental algorithm presented here, is conducted by refining different branches of an 
existing tree via CCD computations. 

This Intersect operation can systematically take advantage of equational constraints. 
The problem of making use of equational constraints in CAD has been studied by many 
researchers (THl [3 [SO] . In Section |6j we provide a detailed discussion on how we solve 
this problem. 

When applied to a polynomial system having finitely many complex solutions, our incre- 
mental CCD algorithm specializes into computing a triangular decomposition, say T), such 
that the zero sets of the output regular chains are disjoint. Moreover, such a decomposi- 
tion has no critical pairs in the sense of the equiprojectable decomposition algorithm of |19j . 
This implies that only the "Merge" part of the "Split & Merge" algorithm of is required 
for turning D into an equiprojectable decomposition (which is a canonical representation of 
the input variety, once the variable order is fixed) . Consequently, one could hope extending 
the notion of equiprojectable decomposition (and related algorithms) to positive dimension 
by means of our incremental CCD algorithm. This perspective can be seen as an indirect 
application of CAD to triangular decomposition. 

As we shall review in Section [2j a CCD is encoded by a tree data-structure. Then each 
path of this tree is a simple system in the sense of [Ml [SZ] ■ So the work presented here can 
also be used to compute a Thomas decomposition of a polynomial system [371 [5] • Moreover, 
the decomposition we compute is not only disjoint, but also cylindrically arranged. 

The complexity of our algorithm cannot be better than doubly exponential in the num- 
ber of variables [6] . So the motivation of our work is to suggest possible ways to improve 
the practical applicability of CAD. The benchmark in Section [7] shows that TCAD outper- 
forms Qepcad [MllS] and Mathematica [33] for many well-known examples. The algorithm 
presented in this paper can support QE. We have realized a preliminary implementation of 
an algorithm for doing QE via TCAD. We will report on this work in a future paper. 
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2 Complex cylindrical tree 



Throughout this paper, we consider a field k of characteristic zero and denote by K the 
algebraic closure of k. Let k[x] be the polynomial ring over the field k with ordered variables 
X = < • • • < Xn- Let p S k[x] be a non-constant polynomial and a; £ x be a variable. 
We denote by deg{p,x) and lc{p,x) the degree and the leading coefficient of p w.r.t. x. 
The greatest variable appearing in p is called the main variable, denoted by mvar(p). The 
leading coefficient, the degree, the reductum of p w.r.t. mvar(p) are called the initial, the 
main degree, the tail of p; they are denoted by init(p), mdeg(p), tail(p) respectively. The 
integer k such that xj. = mvar(p) is called the level of the polynomial p. We denote by 
der(p) the derivative of p w.r.t. mvar(p). The notions presented below were introduced 
in fT¥| and they are illustrated at the beginning of Section [sj 

Separation. Let C be a subset of K"~^ and P C k[a;i, . . . , Xn-i,Xn] be a finite set of level 
n polynomials. We say that P separates above C if for each a G C: 

• for each p & P, the polynomial init(p) does not vanish at a, 

• the polynomials p{a,Xn) € K[x„], for all p d P, are squarefree and coprime. 
Note that this definition allows C to be a semi-algebraic set, see Theorem [3j 

Cylindrical decomposition. By induction on n, we define the notion of a cylindrical 
decomposition of K" together with that of the tree associated with a cylindrical decompo- 
sition of K". For n = 1, a cylindrical decomposition of K is a finite collection of sets 
2? = {Di, . . . , Dr+i}, where either r = and Di = K, or r > and there exists r non- 
constant coprime squarefree polynomials pi, . . . ,pr of k[xi] such that for 1 < z < r we have 
Di = {xi € K I Pi{xi) — 0}, and Dr+i = {xi G K | pi{xi) ■ ■ -prixi) ^ 0}. Note that the 
Di's, for all 1 < « < r -I- 1, form a partition of K. The tree associated with I? is a rooted tree 
whose nodes, other than the root, are Di, . . . , D^, D^+i which all are leaves and children 
of the root. Now let n > 1, and let V — {Di, . . . , Dg} be any cylindrical decomposition 
of K"~^. For each Di, let be a non-negative integer and let {pi^i, ■ ■ ■ ,Pi.ri} be a set of 
polynomials which separates above Di. If = 0, set Di^i = x K. If > 0, set 



for 1 < j < r.i and set = G K" | a G A and [YYj^iPi,jia,x„)] ^ 



The collection V = {Dij | 1 < i < s, 1 < j < + 1} is called a cylindrical decomposition 
of K". The sets Dij are called the cells of V. If T' is the tree associated with V' then the 
tree T associated with V is defined as follows. For each 1 < i < s, the set Di is a leaf in T' 
which has all Dij's for children in T; thus the Aj's are the leaves of T. 

Note that each node of T is either associated with no constraints, or associated with 
a polynomial constraint, which itself is either an equation or an inequation. Note also that, 
if the level of the polynomial defining the constraint at A'' is i, then i is the length of a path 
from N to the root. Moreover, the polynomial constraints along a path from the root to 
a leaf form a polynomial system called a cylindrical system o/k[xi, . . . ,Xn] induced by T. 
Let S be such a cylindrical system. We denote by Z{S) the zero set of S. Therefore, each 
cell of V is the zero set of a cylindrical system induced by T. 

Let r be a sub-tree of T such that the root of F is that of T. Then, we call F a 
cylindrical tree of k[xi, . . . , a;„] induced by T. This cylindrical tree F is said partial if it 
admits a non-leaf node N such that the zero set of the constraint of N is not equal to the 
union of the zero sets of the constraints of the children of N . If F is not partial, then it is 
called complete. 



A J = {{a,Xn) G K" I a G A and a;„) = 0}, 
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In the algorithms of Section|4j the cyhndrical tree is an essential data structure. SectionjS] 
discusses the main properties and operations on this data structure. 

Let F = {/i, . . . , fs} be a finite set of polynomials of k[a;i < • • • < A cyhndrical 
decomposition V of K" is called F -invariant if for any given cell D oi T> and any given 
polynomial f £ F, either / vanishes at all points of Z) or / vanishes at no points of D. 

Example 1. Let F :~ {y^ +x,y'^ +2/}- F-invariant cylindrical decomposition of is 
illustrated by Figure^ 




y'- + y^Q y+l = Q ]/' - y ^ Q y'- + y = 

Figure 1: an F {y^ + + y} invariant complex cylindrical tree 



We observe that every cylindrical system induced by a cylindrical tree is a simple system, 
as defined by Wang in [37]. This notion was first introduced by Thomas in 1937 [35]. Simple 
systems have many nice properties. For example, if [A, B] is a simple system, then the pair 
[A, HpesP] ^ squarefree regular system, as defined by Wang in [371 [35] . 

Let r be a cylindrical system of k[x] and let p be a polynomial of k[x]. We say that p 
is invertible modulo T if for any a S Z{T), we have p{a) ^ 0. We say that p is zero modulo 
r if for any a € ■^(F), we have pia) = 0. We say that p is sign invariant above F if p is 
either zero or invertible modulo F. Let q be another polynomial of k[x]. We say that p = q 
modulo F if Z{T) n Z{p) = Z{T) f] Z{q). 

Greatest common divisor (GCD). Let p and / be two level n polynomials in k[x]. Let F 

be a cylindrical system of k[xi, . . . , Xn-i]- For any u G K"~^ of Z{r), assume at least one 
of lc(p, x„)(m) and lc(/, is not zero. A polynomial g e k[x] is called a GCD of p and 

f modulo V if for any u € K"-i of Z(F), 

• g{u) is a GCD of p{u) and f{u) in K[a;„], and 

• we have \c{g,Xn){u) ^ 0. 

Let dp = deg(p, Xn), df — deg(/, x„). Recall that we assume dp,df > 1. Let A = 
min(cip,d/). Let F be a cylindrical system of k[xi, . . . , x„_i]. Let 5*0, . . . , S'a-i be the 
subresultant polynomials |31l I22j of p and / w.r.t. x„. Let Si = coeff(5i, xjj be the princi- 
ple subresultant coefficient of S',, for < i < A — 1. If dp > d/, we define S\ — /, S'a+i = p, 
sx = init(/) and sa+i = init(p). If dp < df, we define S\ = p, S\+i — f, s\ — init(p) and 
SA+i = init(/). 



4 



Theorem 1. Let j be an integer, with 1 < j < X + I, such that Sj is invertible modulo T 
and such that for any < i < j , we have s,; = modulo T. Then Sj is a GCD of p and f 
modulo r. 

Proof. It can be easily proved by the specialization property of subresultant chains. In 
particular, it is a direct corollary of Theorem 5 in [13] . □ 



3 Data structure for cylindrical decomposition 

In this section, we describe the data-structures that are used by the algorithms presented 
in this paper for computing cylindrical decompositions. To understand the motivation of 
our algorithm design, let us consider a simple example with n — 2 variables. Let a, b 
be two coprime squarefree non-constant univariate polynomials in fc[a::i]. Observe that 
L := k[xi]/ (ab) is a direct product of fields. Let also c, d be two bivariate polynomials of 
k[xi,X2], such that deg(c, a;2) > 0, deg{d,X2) > 0, and lc(c, 0:2) = lc{d,X2) = 1 hold and 
such that c, d are coprime and squarefree univariate as polynomials of iv[x2]. Therefore the 
following four polynomial systems are simple systems 

a{xi)b{xi) = ( a{xi)b{xi) ^ ( a{xi)b{xi) ^ S -r \h( ^ \ ^ (\ 

C(X1,X2)=0 '\ d{x^,X2)^0 ' \ c{xi,X2)d{xi,X2)^Q ' ^ ^^^^ ^'^ 

that we denote respectively by 5*1 , S'2 , /Ss , S'4 . It is easy to check that the zero sets Z{Si), 
Z{S2), Z{S3), Z{S4) are the cells of a cylindrical decomposition V of K^. 

Let / G k[xi] be another univariate polynomial. Assume that one has to refine V 
into a cylindrical decomposition of which is required to be {/}-invariant. That is, one 
has to test whether / is invertible or zero modulo each of the systems 6*1 , ^2 , 5*3 , S'4 , and 
further decompose when appropriate. Assume that the polynomial a divides / whereas 6, / 
are coprime. Assume also that the system Si is processed first in time. By computing 
gcd(/, a6), which yields a, one splits Si into the following two sub-systems that we denote 
by Si^i and Si^2- 

= r b{xi) = 

c{xi,X2)^0 ' \ c(.Xi,X2)=0. 

Assume that S2 is processed next. By computing gcd(/, ab) (again) one splits S2 into the 
following two sub-systems that we denote by S'2,1 and S'2,2- 

aixi)^0 / 6(xi)=0 



d{xi,X2)=0 ' \d{xi,X2)^0- 

Consequently, in the course of the creation of 5i.2, >5'2,i and S2^2y the same polynomial 
GCD and the same field extensions (namely k[xi]/ (a) and k[xi]/{b)) were computed twice. 
This duplication of calculation and data is a common phenomenon and a performance 
bottleneck in most algorithms for decomposing polynomial systems. 

Mathematically, each constructible set should not be represented more than once in a 
computer program. To implement this idea, all constructible sets manipulated during the 
execution of a given computer program should be seen as part of the same universe, say K". 
Moreover, the subroutines of this program should have the same view on the universe, which 
is then a shared data- structure, such that whenever a subroutine modifies the universe all 
subroutines have immediate access to the modified universe. Satisfying these requirements is 
a well-known challenge in computer science, an instance of which is the question of memory 



5 



consistency for shared-memory parallel computer architectures, such as multicores. With 
our above example, even if we do not intend to run computations concurrently, we are 
concerned with the practical efficiency and ease-of-use of the mechanisms that maintain 
up-to-date all views on the universe. 

Recall that a cylindrical decomposition can be identified to a tree where each node is a 
constructible set of K" given by either an equation constraint, or an inequation constraint, 
or no constraints at all. In this latter case, the corresponding constructible set is the whole 
space. All algorithms in Sectionj4]work on a given cylindrical decomposition D encoded by 
a tree T (as defined in Section [2^7 That is, the tree T is regarded as the universe. 

We assume that there is a procedure for updating the tree T, which, given a "node- 
to-be-replaced" N and its "replacing nodes" A^i, . . . , Ne, is called split(A^; A^i, . . . , N^) and 
works as follows: 

1. for i — 1, . . . ,e, for each child C of deeply copy (thus creating new nodes) the 
sub- tree rooted at C and make that copy of C a child of Ni, 

2. update the parent of such that Ni, . . . , are new children of the the parent of A^, 

3. remove the entire sub-tree rooted at A^ from the universe, including A^. 

We assume that all updates are performed sequentially (thus using mutual exclusion mecha- 
nism in case of concurrent execution of the algorithms of Section |4]) such that no data-races 
can occur. 

We also assume that each node A^ (whether it is a node in the present or has been re- 
moved from the universe) has a unique key, called key(A^), and a data field, called value(A^), 
storing various information including: 

• a time stamp PAST or present, 

• if PAST, the list of its replacing nodes (as specified with the split procedure) and the 
list of its children at the time it was replaced, 

• if PRESENT, the list of its children and a pointer to the parent. 

All nodes are stored in a dictionary H which can be accessed by all subroutines. Mod- 
ifying the universe means updating H using the split procedure. Since all our algorithms 
stated in Section [4] are sequential, no synchronization issue has to be considered. The 
mechanism described above allows us to achieve our goals. 

4 Constructing a cylindrical tree incrementally 

In this section, we present an incremental algorithm for computing a cylindrical tree, as 
defined in Section [2] We start by commenting on the style of the pseudo-code. Secondly, 
we present the specifications of the algorithm and related subroutines. Thirdly, we state all 
the algorithms in pseudo-code style. Finally, proof sketches of the algorithms are provided 
at the end of this section. 

Following the principles introduced in Section [3j our procedures operate on a "universe" 
(which is a cylindrical tree T) that they modify when needed. These modifications are of 
two types: 

• splitting a node, 

• attaching information to a node. 

In addition to the attributes described in Section |3j a node has attributes corresponding 
to the results of operations like Squarefree, Gcd, Intersect. In other words, our procedures 
do not return values; instead they store their results in the nodes of the universe. This 
technique greatly simplifies pseudo-code. 
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Since attributes of nodes are intensively used in our pseudo-code, we use the standard 
"dot" notation of object oriented programming languages. In addition, since a node can 
have many attributes, we make the following convention. Suppose that a node V is split 
into two nodes Vi and V2. Some attributes are likely to have different values in Vi (resp. V2) 
and V. But most of them will often have the same values in both nodes. Therefore, after 
setting up the values of the attributes that differ, we simply write Vi. others := y. others to 
define the attributes of Vi whose values are unchanged w.r.t. V. 

Several procedures iterate through all the paths of the universe T. By path, we mean 
a path (in the sense of graph theory) from the root of T to a leaf of T. The current 
path is often denoted by F or C. Recall from Section [2] that a path in T corresponds to 
a simple system, say S. Computing modulo S may split S and thus modify the universe 
automatically, that is, in a transparent manner in the pseudo-code. However, splitting S 
also changes the current path. For clarity, we explicitly invoke a function called UpdatePath, 
which updates its first argument (namely the current path) from the universe. 

In order to iterate through all the paths of the universe T, we use a function NextPathToDo. 
This command is a generator or an iterator in the sense of the theory of programming lan- 
guages. That is, it views T as a stream of paths and returns the next path-to-be-visited, if 
any. Thanks to the fact that the universe is always up-to-date, the function NextPathToDo 
is able to return the next path-to-be-visited in the current state of the universe. 

A frequently used operation on the universe and its paths is ExtractProjection, see for 
instance Algorithm [6j When applied to the universe T and an integer k (for < fc < n, 
where n is the length of a path from the root of T to a leaf of T) ExtractProjection returns 
a "handle" on the universe "truncated" at level k, that is, the universe where all nodes of 
level higher than k are ignored (thus viewing the level k nodes as leaves). When applied to 
path, ExtractProjection has a similar output. 

We often say that a function (see for instance Algorithm [s]) returns a refined cylindrical 
decomposition. This is another way of saying that the universe is updated to a new state 
corresponding to a cylindrical decomposition refining (in the sense of a partition of a set 
refining another partition of the same set) the cylindrical decomposition of the previous 
state. 

After these preliminary remarks on the pseudo-code, we present the specifications of the 
algorithm and related subroutines. 

The top level algorithm for computing a cylindrical tree is described by Algorithm [4] 
It takes a set F of non-constant polynomials in k[xi < • • • < a;„] as input and returns an 
_F-invariant cylindrical decomposition of K". This algorithm relies on a core operation, 
called Intersect, which computes a cylindrical decomposition in an incremental manner. 

The Intersect operation is described by Algorithm [Sj It takes a cylindrical tree T and 
a polynomial p of k[2;i < • • • < Xn] as input. It refines tree T such that p is sign invariant 
above each path of the refined tree T. This operation is achieved by refining each path of 
T with IntersectPath. 

The IntersectPath operation is described by Algorithm |6] It takes a polynomial p, a 
cylindrical tree T and a path F of T in k[a;i < • • • < x„] as input. It refines F and updates 
the tree T accordingly such that p is sign invariant above each path derived from F in the 
updated tree T. This operation finds the node N in F whose level is the same as that of 
p. Let Fjv be the sub-path of F from N to the root of T. The IntersectPath operation then 
calls the routine IntersectlVlain so as to refine F^r into a tree T/v such that p becomes sign 
invariant w.r.t. Tjy. 

The routine IntersectMain is described by Algorithm [7] It takes a cylindrical tree T, a 
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path r of r, and a polynomial of the same level as the leaves of T in k[a;i < • • • < Xn] as 
input. It refines T and updates the tree T accordingly such that p becomes sign invariant 
above each path derived from F in the updated tree. 

The routine IntersectMain works in the following way. It first splits F such that above 
the projection Cn-i of each new branch C of F in K"~^, the number of distinct roots 
of p w.r.t. Xn is invariant. This is achieved by the operation Squarefree, described by 
Algorithm [8] The squarefree part of p above a branch C is denoted by sp. If p has no 
roots or is identically zero above Cn-i, the sign of p above C is determined immediately. 
Otherwise, a case discussion is made according to the structure of the leaf node V of C. 
If V has no constraints associated to it, then V is simply split into two new nodes sp — 
and sp 7^ 0. Assume now that V has a constraint, which can be either of the form / = or 
of the form / ^ 0, where / is a level n polynomial squarefree modulo C„_i. This case is 
handled by computing the GCD g of sp and / modulo C„_i. The node V then splits based 
on the GCD g and the co- factors of sp and /. 

The GCD is computed by the operation Gcd, described by Algorithm [9| and [TO) The 



CO- factors are computed by Algorithm 11 The Squarefree and Gcd operations rely on the op- 
eration MakeLeadingCoefficientlnvertible, described by Algorithm [T2| This latter operation 
takes as input a polynomial p ofh[xi < ■ ■ ■ < Xn] , a cylindrical tree T of k[a;i < • • • < 
and a path F of T. Then, it refines F and updates T accordingly such that above each path 
C of T derived from F, the polynomial p is either zero or its leading coefficient is invertible. 

All the algorithms also rely on the following three operations which perform manipula- 
tions and traversal of the tree data structure. For these three operations, only specifications 
are provided below while their algorithms are explained in Section [3j 

Algorithm 1: UpdatePath(F, T) 

- Input: A cylindrical tree T. A path F in some past state of T. 

- Output: A subtree ST in present state of T. ST is derived from F according to the 
historical data of T. 



Algorithm 2: ExtractProjection(T, fc) 

- Input: A cylindrical tree T of k[a::i < • • • < x„]. An integer k, < k < n. 

- Output: A cylindrical tree Tk in k[.Ti < • • • < Xk] such that Tk is the projection of T 
in k[xi < ■ ■ ■ < Xk]- 



Algorithm 3: NextPathToDo„(T) 

- Input: A cylindrical tree T in k[a;i < • • • < x„]. 

- Output: For a fixed traversal order of a tree, return the first "ToDo" path F of T. 



Theorem 2. For a set of polynomials in k[xi, . . . , a;„], Algorithm^ computes an F- 
invariant cylindrical decomposition of K" . 

Proof. Firstly, we prove the termination. The basic mutual calling graph of its subroutines 
are: 

lntersectMain„ SquarefreCn — > lntersectMain„_i 
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Algorithm 4: CylindricalDecompose(F) 



Input: F is a set of non-constant polynomials in k[a:i < • • • < a;„]. 
Output: An F-invariant cylindrical decomposition of K". 
1 begin 

2 

3 

4 
5 



create a tree T with only one vertex Vq: the root of T; 
for i from 1 to n do 

create a vertex Vf, Vi. signs := 0; Vi. formula := "any x"; 

Vi-^i.child := Vi; 

for p G F do 
|_ lntersect„(p,r); 

return T; 
9 end 



Algorithm 5: lntersect„(p, T) 

Input: A cylindrical tree T of k[xi < • • • < a;„]. A non-constant polynomial p of 

k[a;i < • • • < Xn]- 

Output: A refined cylindrical decomposition such that p is sign invariant above 
each path of T. 

1 while r := NextPathToDo„(T) ^ do 

2 [_ lntersectPath„(p,r,T); 



Algorithm 6: lntersectPath„(p, T, T) 



Input: A cylindrical tree T of k[a;i < • • • < a;„]. A path F of T. A polynomial p of 

k[.Ti < ■■ ■ < Xn]. 

Output: A refined cylindrical decomposition T such that p is sign invariant above 
each path derived from F. 

1 begin 

if p e k then 
I return; 
else 

k := level (p); 
if fc = n then 
I lntersectMain„(p,r,T); 
else 

Tk := ExtractProjection(T, A;);Ffe := ExtractProjection(F, fc); 
lntersectMainfc(p, F^, T^); 
UpdatePath(r,T); 
for each leaf V of F do 

1^ Let Lk be the ancestor of V of level k; V.signs\p] := Li..signs\p] ; 



14 end 
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Algorithm 7: lntersectMain„(p, T, T) 



Input: A cylindrical tree T of k[xi < • • • < a;„]. A path F of T. A polynomial p of 

level n in k[a;i < • • • < x„]. 
Output: A refined cylindrical decomposition T such that p is sign invariant above 
each path derived from T. 

begin 

T„_i := ExtractProjection(r, n — 1); r„_i := ExtractProjection(r, n — 1); 
Squarefree„(p, r„_i, T„_i); 
UpdatePath(r,r); 

while C := NextPathToDo„(r) ^ do 

V := C.leaf; C„_i := ExtractProjection(C, n — 1); 

sp := Cn-i-leaf .Square free\p]; 
if = then 
I V.signs\p] := 0; 
else if sp = 1 then 
I V.signs\p] := 1; 

else if V. formula is "any Xn " then 

split V into two new vertices Vi and V2; 
Vi. formula := sp = 0; Vi. signs := V.signs; Vi.signs\p] := 0; 
V2-formula := sp ^ 0; V2. signs := V.signs; V2.signs\p] := 1; 
Vi. others := V. others; V2.others := V.others; 
Cn-i -leaf .children := Vi,V2; 



1 
2 
3 

4 
5 

6 
7 
8 
9 
10 
11 
12 
13 
14 
15 
16 
17 
18 

19 
20 
21 
22 
23 
24 
25 
26 
27 
28 
29 
30 
31 
32 
33 
34 

35 
36 
37 
38 
39 
40 
41 
42 
43 
44 



else 



// V. formula is of the form f = or f^O 
Gcdn{sp, /,C„_i,r„_i); 
UpdatePath(C,T); 
for each leafV of C do 
let L be the parent of V; 
cp, g, cf := CoFactor(sp, L.Gcd[sp, /],/); 
if V.form,ula is of the form / = then 
if g = 1 then 
I V.signs\p\ := 1; 
else if cf = 1 then 
I V.signs[p] := 0; 
else 

split V into two new vertices Vi and V2; 
Vi. formula .= g = 0;Vi. signs := V.signs; Vi.signs\p] := 0; 
V2. formula := cf = 0;V2. signs := V.signs; V2.signs\p] := 1; 
Vi. others := V.others; V2.others := V.others; 
L.children := Vi, V2; 



else 



if cp = 1 then 

I V.signs\p] := 1; 
else 

split V into two new vertices Vi and V2; 
Vi.form,ula := cp = 0; Vi. signs := V.signs; Vi.signs\p] := 

V2. formula := (/ * cp) 7^ 0; 
V2 .signs := V.signs; V2.signs[p] := 1; 
Vi.others V.others; V2-others := V.others; 
L.children := Vi, V-^q 



0; 



45 end 



Algorithm 8: Squarefree„(p, F, T) 



Input: A cylindrical tree T of k[a;i < • • • < Xn-i]- A path F of T. A polynomial p 
of level n. 

Output: A refined cylindrical tree T of k[.xi < • • • < .t„-i]. Above each path C of 
T derived from F, there is a dictionary C. leaf .Square free. Let 
p* := C.leaf. Square free\p]. We have: 

• p = p* modulo C . 

• If p* is of level n, then both init(p*) and discrim(y)*) are invertible modulo C. 

• If p* is of level less than n, then p* is either or 1. 

1 begin 

2 if n = 1 then 

3 let r be the root of T; r. Square free\p] := SquarefreePart(p); 

4 return 

5 MakeLeadingCoefficientlnvertible„(p,p, F,T); 

6 while C := NextPathToDo„_i(F) ^ do 

7 / := C.leaf .InvertLc[p]; 

8 if level(/) <n or deg(/, = 1 then 

9 I C .lea f .Square free[p\ := / 

10 else 

11 Gcd„(/,der(/),C,r); 

12 for each leaf L of C do 

13 g:= L.Gcd[/, der(/)]; 

14 if g = 1 then 

15 I L.Squarefree[p] := f 

16 else 

17 1^ L. Square free\p\ :=pquo(/, 5) 



18 end 



Algorithm 9: Gcd„(p, /, F, T) 



Input: A cylindrical tree T of k[xi < • • • < Xn-i]- A polynomial 

p € k[a;i < • • • < Xn] of level n. A path F of T. A polynomial / of level n 
such that init(/) is invertible modulo F. 
Output: A refined cylindrical tree T. Above each path C of T derived from F, 

there is a dictionary C.leaf .Gcd such that C.leaf .Gcd[p, f] is a GCD of p 
and / modulo C. 

1 begin 

2 let S be the subresultant chain of p and /; 

3 if mdeg(j?) > mdeg(/) then 

4 \ d:= mdeg(/) 

5 else 

1^ d := mdeg(p) + 1 

return Gcd„(p, /, S, d, 0, F, T); 
8 end 
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Algorithm 10: Gcdn{pJ,S,d,i,T,T) 



Input: 

• A polynomial p G k[a;i < • • • < a;„] of level n. 

• A polynomial / of level n such that lc(/) is invertible modulo F. 

• The subresultant chain S of p and / w.r.t. Xn- 

• A non-negative integer d (as defined in the pseudo-code of Algorithm |9]) and such 
that the principle subresultant coefficient is invertible modulo T. 

• A non-negative integer i such that < i < d and the principle subresultant 
coefficient Sj is zero modulo F, for all < j < z. 

• A path F of T. 

• A cylindrical tree T of k[a;i < • • • < a;„_i]. 

Output: A refined cylindrical tree T. Above each path C of T derived from F, 

there is a dictionary C.leaf.Gcd such that C.leaf.Gcd[p, f] is a GCD of p 
and / modulo C. 

begin 

i{ i — d then 

T.leaf.Gcd[p, /] := Si; 
return; 

lntersectPath„_i(si, F, T); 
while C := NextPathToDo„_i(F) 7^ do 
if C.leaf.signs[si] — 1 then 
if i = then 
I C.leaf.Gcd[p, f] := 1 
else 

I C.leaf.Gcd[p, f] := 

else 

I Gcd„(p,/,^,d,i + l,C,T) 



7 

8 
9 
10 
11 

12 
13 



14 end 
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Algorithm 11: CoFactor(p, g, /) 



Input: Two polynomials p and / of level n in k[a;i < • • • < A polynomial g 

which is cither 1 or of level n in k[a;i < • • • < a;„]. 
Output: As described by the algorithm. 
1 begin 



it g = 1 then 

I cp:=p; gg := 1; cf := /; 
else if mdeg(gf) = mdcg(/) then 
.9.9 := ./; 

if mdcg(g) = mdeg(p) then 

I c.f 1; cp := 1; 
else 

|_ cf := 1; cp := pquo{p,gg) 

else if mdeg(g) = mdeg(p) then 
I 99 ■= P; cf := pquo(/, gg); cp := 1; 
else 

|_ cp := pqno{p,g); cf := pquo{f,g); gg := g; 
return cp, gg, cf; 



15 end 



Algorithm 12: MakeLeadingCoefficientlnvertible„(p,p, F, T) 

Input: A polynomial p of k[.Ti < • • • < x„]. A polynomial p of k[a:i < • • • < Xn] 
such that p = p modulo T. A cylindrical tree T of k[xi < ■ ■ < Xn-i]- A 
path r of T. 

Output: A refined cylindrical tree T of k[.xi < • • • < x„_i]. Above each path C of 
T derived from T, there is a dictionary C.leaf.InvertLc. Let p* be the 
polynomial C. leaf .Invert Lc\p]. Then, we have: 

• p = p* modulo C. 

• If p* is of level n, then init(p*) is invertible modulo the path C. 

• If p* is of level less than n, then p* is either or 1. 
begin 

lntersectPath„_i (lc(p, F, T); 
while C := NextPathToDo„_i(F) ^ do 

if C.leaf.signs[lc{p,Xn)] = 1 then 
if level(p) < n then 
I C.leaf.InvertLc\p] := 1 
else 

|_ C.leaf.InvertLc\p] := p 

else 

if level(p) < n then 
I C.leaf.InvertLc\p\ := 
else 



1 

2 
3 
4 
5 
6 
7 
8 

9 
10 
11 
12 
13 



|_ MakeLeadingCoefficientlnvertible„(p, tail(p), C, T) 



14 end 
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and 



lntersectMain„ Gcdp lntersectMain„_i 



So the termination is easily proved by induction. The correctness follows from the specifi- 
cation of its subroutines and Theorem [ij □ 

Example 2. In this example, we illustrate the operation IntersectPath. Let F :— {y^ + 

XjU^ + y}- The incremental algorithm first computes an + x sign invariant complex 
cylindrical tree, which is described by the following tree T . 

„2 ^ ^ _ 



T := 



X 







y = Q 



- X : 

x^O 



x^Q 



a; = 



x = Q 
x^O 



Let r be the path {x — 0,y ^ 0} ofT. Calling lntersectPath(2/^ + y, T, T) will update T into 
the following tree. 



x^Q 



X 7^ 




y2 + 2: = 



y 



x = Q 
x^Q 



5 Building a CAD tree from a complex cylindrical tree 

In this section, we review briefly how to compute a CAD of M" from a cylindrical decom- 
position of C". The reader may refer to [T^ for more details. Recall that n > 1 holds. 
We denote by 7r„_i the standard projection from M" to W^~^ that maps {xi, . . . , Xn-i,Xn) 

onto {xi,...,Xn-l)- 

Stack over a connected semi-algebraic set. Let 5 be a connected semi- algebraic subset 
of M."-^. The cylinder over S in M" is defined as Zm(S') := S x R. Let Bi < ■ ■ ■ < Os 
be continuous semi-algebraic functions defined on S. The intersection of the graph of Oi 
with Z^{S) is called the di-section of Zr(5). The set of points between two consecutive 
sections of Zr(S') is a connected semi-algebraic subset of M", called a sector of Zr(S'). All 
the sections and sectors of Z^{S) form a disjoint decomposition of Z^{S), called a stack 
over S. 

Cylindrical algebraic decomposition. A finite partition T) of R" is called a cylindrical 
algebraic decomposition (CAD) of M" if one of the following properties holds. 

• Either n = 1 and I? is a stack over MP . 

• Or the set of {iTn-iiD) | D e P} is a CAD of W"-^ and each I? e 2? is a section or 
sector of the stack over 7r„_i(Z?). 

When this holds, the elements of T) are called cells. 

Sign invariance and delineability. Let p be a polynomial of M[a:i, . . . ,a;„], and let S' be a 

subset of M". The polynomial p is called sign invariant on S if the sign of does not 
change when a ranges over S. Let F C M.[xi, . . . , a;„] be a finite polynomial set. We say S 
is i^-invariant if each p ^ F in invariant on 5. A cylindrical algebraic decomposition T) is 
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i^-invariant if F is invariant on each cell D (z V. Let p be a polynomial of M[a;i, . . . 
and let S* be a connected semi-algebraic set of M""^. We say that p is delineable on S if the 
real zeros of p define continuous semi- algebraic functions 9i, . . . ,9^ such that, for all a € S* 
we have Oi{a) < ■ ■ ■ < Os(a). In other words, p is delineable on S if its real zeros naturally 
determine a stack over S. We recall the following Theorem introduced in [Tl]. 

Theorem 3. Let P — {pi, . . . ,pr} be a finite set of polynomials in R[xi < ■ ■ ■ < x„] of 
level n. Let S be a connected semi-algebraic subset o/M"~^. If P separates above S, then 
each Pi is delineable on S . Moreover, the product of the pi, . . . ,pr is also delineable on S . 

Let F be a finite set of polynomials in Q[a;i < • • • < a;„]. Let CT be an F-invariant 
complete cylindrical tree of C" . Applying Theorem [3] to polynomials in CT, we can derive 
an F-invariant cylindrical algebraic decomposition of M" by induction on n. A procedure 
MakeSemiAlgebraic, was introduced in [Hj to derive a CAD from a CT via real root isolation 
of zcro-dimensional regular chains. 



Example 3. Let F -.^ {y^ + x}. 

described by the following tree. 



An F -invariant cylindrical algebraic decomposition is 



T := 



a; < 



X = 




x\Ay< yia^ 



y"^ -\~ X > 
y^ -\-x^0 
y"^ -\- X > 



y'^ -\- X > 

y 
y 



+ x^0 
2 < 



y^ + x ^0 
y'^ + X > 



X > for any y : y^ + x > 



6 Making use of equational constraints and other opti- 
mizations 

In this section, we discuss several possible optimizations to algorithms presented in Sec- 
tion m 

Firstly, we discuss how to compute a CAD dedicated to a semi-algebraic system, which 
provides a systematic solution for making use of equational constraints when computing 
CADs. The motivation for making use of equational constraints comes from quantifier 
elimination. Let 

PF := {Qk+lXk+l ■ ■ ■ QnXn)FF{xi, . . . , x„), 

be a prenex formula, where FF is a DNF formula. To perform QE by CAD, the first 
computation step is to collect all the polynomials appearing in FF as a polynomial set F 
and compute an _F- invariant CAD of M" . This process of computing an F-invariant CAD 
exhausts all possible sign combinations of F, including those which do not appear in FF, and 
thus often computes much more than needed for solving the input QE problem. Different 
techniques in the literature have been proposed for taking advantage of the structure of 
the input problem. These methods include partial CAD |17j for lazy lifting, simplified 
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projection operator for handling pure strict inequalities (23l33|, smaller projection sets for 
making use of equational constraints [HJ [3 1301 • 

To make the discussion clear, we first quote a paragraph of [7]. "The idea is as follows: 
if an input formula includes the constraint / = 0, then decompose W into regions in which 
/ has invariant sign, and then refine the decomposition so that the other polynomials have 
invariant sign in those cells in which / = 0. The signs of the other polynomials in cells in 
which / 7^ are, after all, irrelevant. Additionally, the method of equational constraints 
seeks to deduce and use constraints that are not explicit in the input formula, but rather 
arise as consequences of two or more explicit constraints (e.g. if / = and g = are explicit 
constraints, then res(/, g) = is also a constraint.)" 

This idea, of course, is attractive. Much progress on it has also been made. However, 
the reason why it is a generally hard problem for CAD is that the framework of PCAD 
does not have much flexibility to allow propagation of equational constraints. In the world 
of PCAD, one always tries to obtain a generic projection operator and then applies the 
same projection operator recursively. To obtain a generic projection operator for handling 
equational constraints is hard because many problems inherently require different projection 
operators during projection. Therefore case discussion is important. 

In fact, case discussion is very common in algorithms for computing triangular decompo- 
sitions. For such algorithms, equational constraints are natural input of these algorithms. 
The two keys ideas "splitting only above / = 0" and "if / = and 5 = are explicit 
constraints, then res(/, 5) = is also a constraint" have already been systematically taken 
care of in the Intersect operation of the authors' paper for computing triangular decompo- 
sitions [l2] . 

Next we explain how to modify algorithms presented in Section [4] to automatically 
implement these ideas. 

Suppose now that the input of Algorithm CylindricalDecompose is a system of equations 
or inequations, this algorithm will then compute a partial cylindrical tree such that its zero 
set is exactly the zero set of input system. This can be simply achieved by passing an 
equation or inequation to the function Intersect. W.l.o.g., let us assume that an equation 
p = is passed as an argument of Intersect. Then for this function and all its called 
subroutines, we will cut the computation branches above which p is known to be nonzero 
and never proceed with computation branches above which p cannot be zero. For example, 
we will not create a new vertex at step 15, 32, 42 in Algorithm IntersectMain. We will delete 
the vertex V at step 11, 26, 37 since p is nonzero on V. 

The first important optimization in IntersectMain which can be implemented is to avoid 
Squarefree computation at step 3 if T.leaf is an equational constraint. This idea is quite 
close to "splitting only above / = 0" . Another important optimization can be done at step 
19 of IntersectMain. Assume that V. formula is an equational constraint / = 0, then when 
Gcd is called, in step 5 of Algorithm [lO) we can do as follows. If i = 0, then Si is the 
resultant of p and /. Thus we should pass = to the IntersectPath operation in order 
to avoid useless computations on the branch Si ^ 0. This addresses the idea "if f — 
and g = are explicit constraints, then res(/, 5) = is also a constraint." Moreover, these 
optimizations are systematically performed during the whole computation. 

Next we briefly mention several other important optimizations. Let ^ be a leaf of a 
path r of a cylindrical tree. Assume that V. formula is of the form / 7^ or of the form 
/ = 0. We can safely replace / by its primitive part since lc(/) is invertible modulo r„_i. 
Replacing / by its irreducible factors over Q is often a more efficient choice. Last but not 
least, recall that a path F in the cylindrical tree is a simple system. Writing F as two parts 
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r := [T,H], where T is a set of equations and H is a, set of inequations. We know that T 
is a regular chain and F is a squarefree regular system. Thus the Zariski closure of F is the 
variety of the saturated ideal of T. We can call the pseudo division operation preni(p, T) 
or preni(/, T) to test whether p or / is zero modulo F. And sometimes replacing p by 
prem(p, T) and / by prem(/, T) also ease the computations. 

Example 4. Let F :— {y^ + x — 0,y^ + y — 0} be a system of equations. Taking F as 
input, Algorithm CylindricalDecompose generates the following partial cylindrical tree T of 
such that the zero set of F is exactly the union of the zero sets of the paths in T. 




x = Q a; + l =0 




y=0 y+l=0 
Figure 2: A partial cylindrical tree T adapted to F 



7 Benchmark 

In this section, we report on the experimental results of a preliminary implementation in 
the RegularChains library of MAPLE of the algorithms of Sections |4] and [5j 

The examples in Table [l] and Table [2] are from papers on polynomial system solving, 
such as [HI |3] and the references therein. All the tests were launched on a machine with 
Intel Core 2 Quad CPU (2.40GHz) and 8.0Gb total memory. The time-out is set as 1 hour. 
In the tables, the symbol > Ih means time-out. 

The Maple functions are launched in Maple 15 with the latest RegularChains library. 
The memory usage is limited to 60% of total memory. The software Qepcad is launched 
with the option -fA^500000000-|-L200000, where the first option specifies the memory to be 
pre-allocated (about 23% of total memory for our machine) and the second option specifies 
the number of prime numbers to be used. 

In Table [l} we report on timings for computing cylindrical decomposition of the com- 
plex space with different algorithms and options. Each input system is a set of polynomials. 
The notation tcd-rec denotes an implementation of the original recursive algorithm in |14j , 
while the notation tcd-inc denotes the incremental algorithm presented in Section |4] Both 
tcd-rec and tcd-inc take a set of polynomials as input. The notation tcd-eqs refers to an 
optimized version of tcd-inc which makes use of equational constraints, as explained in Sec- 
tion |6] With the implementation tcd-eqs, every input polynomial set is regarded as a set 
of equations (equating each input polynomial to zero). As we can see in Table [l] the incre- 
mental algorithm presented in this paper is much more efficient than the original recursive 
algorithm. The timings of tcd-eqs show that the optimizations presented in Section [6] for 
making use of equational constraints are very effective. 
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Tabic 1: Tiniiugs for coiiipiiliug (:\ liii(lrical (kx-oiuixwitioii of xlic C(juii)lcx sijacc 



System 


tcd-rec 


tcd-inc 


tcd-eqs 


System 


tcd-rcc 


tcd-inc 


tcd-eqs 


AlkashiSinus 


3373.966 


14.568 


4.168 


MontesSlO 


> Ih 


> Ih 


2.952 


Alonso 


9.636 


1.404 


0.700 


MontesS12 


> Ih 


> Ih 


7.528 


Arnborg-Lazard-rev 


2759.940 


2419.543 


16.233 


MontcsS15 


> Ih 


> Ih 


77.048 


Barry 


39.346 


1.808 


0.556 


MontesS16 


> Ih 


> Ih 


8.228 


blood-coagulation-2 


235.310 


9.472 


0.808 


MontesS4 


556.390 


102.122 


0.488 


Bronstcin-Wang 


255.427 


35.990 


1.120 


MontcsS5 


1449.810 


119.059 


1.004 


cdc2-cyclin 


> Ih 


68.920 


65.976 


MoiitesS7 


> Ih 


> Ih 


1.060 


circles 


276.389 


2.280 


0.520 


MontesS9 


269.636 


4.212 


0.980 


genLinSyst-3-2 


916.245 


19.537 


1.384 


nql-5-4 


> Ih 


1.056 


0.528 


genLinSyst-3-3 


> Ih 


160.406 


12.408 


r-5 


68.364 


3.232 


0.876 


Gerdt 


> Ih 


> Ih 


1.188 


r-6 


1456.883 


46.458 


1.200 


GonzalezGonzalez 


141.072 


53.451 


0.732 


Raksanyi 


1471.351 


118.227 


1.000 


hereman-2 


> Ih 


40.042 


0.908 


Rose 


> Ih 


51.855 


1.072 


lhlp5 


31.069 


3.984 


0.648 


Wang93 


> Ih 


> Ih 


18.877 


Maclane 


> Ih 


> Ih 


6.420 


YangBaxter Rosso 


54.895 


1.560 


0.844 



Table 2: Timings for computing CAD 



System 


qepcad 


qepcad-eqs 


mathematica-eqs 


tcad 


tcad-eqs 


Alonso 


7.516 


5.284 


0.74 


61.591 


5.776 


Arnborg-Lazard-rev 


> Ih 


> Ih 


0.952 


> Ih 


17.325 


Barry 


Fail 


216.425 


0.032 


8.580 


1.004 


blood-coagulation-2 


> Ih 


> Ih 


> Ih 


985.709 


7.260 


Bronstein-Wang 


> Ih 


> Ih 


26.726 


333.892 


2.564 


cdc2-cyclin 


> Ih 


> Ih 


0.208 


574.127 


503.863 


circles 


21.633 


5.996 


41.211 


> Ih 


40.902 


GonzalezGonzalez 


10.528 


10.412 


0.012 


214.213 


1.136 


lhlp2 


960.756 


5.076 


0.016 


3.124 


0.952 


lhlp5 


10.300 


10.068 


0.016 


35.338 


1.084 


MontesS4 


> Ih 


> Ih 


0.004 


2682.391 


0.888 


MontesS5 


Fail 


Fail 


> Ih 


> Ih 


9.400 


nql-5-4 


93.073 


5.420 


1303.07 


113.675 


1.004 


r-5 


> Ih 


1802.676 


0.016 


1282.928 


1.208 


r-6 


> Ih 


> Ih 


0.024 


> Ih 


1.500 


Rose 


Fail 


> Ih 


> Ih 


606.361 


3.136 


AlkashiSinus 


> Ih 


> Ih 


2.232 


> Ih 


58.775 


genLinSyst-3-2 


Fail 


Fail 


217.062 


3013.764 


6.588 


MontesSlO 


> Ih 


> Ih 


> Ih 


> Ih 


22.797 


MontesS12 


> Ih 


> Ih 


> Ih 


> Ih 


330.996 


MontesS15 


> Ih 


> Ih 


0.004 


> Ih 


395.964 


MontesS7 


> Ih 


> Ih 


245.807 


> Ih 


2.452 


MontesS9 


Fail 


Fail 


> Ih 


110.902 


4.944 


Wang93 


Fail 


Fail 


> Ih 


> Ih 


152.673 
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In Table [2] we report on timings for computing CAD with three different computer 
algebra packages: Qepcad, the CylindricalDecomposition command of Mathematica and the 
algorithm presented in Section [4] Each system is a set of polynomials. Two categories 
of experimentation are conducted. The first category is concerned with the timings for 
computing a full CAD of a set of polynomials. For Mathematica, we cannot find any options 
of CylindricalDecomposition for computing a full CAD of a set of polynomials. Therefore for 
this category, only the timings of Qepcad and TCAD are reported. The second category 
is concerned with the timings for computing a CAD of a variety. For this category, the 
timings for Qepcad, Mathematica and TCAD are all reported. 

The notation qepcad denotes computations that Qepcad performs by (1) treating each 
input system as a set of non-strict inequalities and, (2) treating all variables as free vari- 
ables and, (3) executing with the "full-cad" option. The notation tcad corresponds to 
computations that TCAD performs by (1) treating each input system as a set of non-strict 
inequalities and, (2) computing a sign invariant full CAD of polynomials in the input sys- 
tem and, (3) selecting the cells which satisfy those non-strict inequalities. In this way, both 
qepcad and TCAD compute a full CAD of a set of polynomials. 

The notation qepcad-eqs denotes the computations that Qepcad performs by (1) treat- 
ing each input system as a set of equations and, (2) treating all variables as free variables 
and, (3) executing with the default option. The notation mathematica-eqs represents com- 
putations where the CylindricalDecomposition command of Mathematica treats each input 
system as a set of equations. The notation tcad-eqs corresponds to computations where 
TCAD treats each input system as a set of equations. 

From Table [2j we make the following observations. When full CADs are computed, 
within one hour time limit, Qepcad only succeeds on 6 out of 24 examples while TCAD 
succeeds on 14 out of 24 examples. When CADs of varieties are computed, for all the 10 
out of 24 examples that Qepcad can solve within one hour time limit, both Mathematica 
and TCAD succeed with usually less time. For the rest 14 examples, TCAD solves all of 
them while Mathematica only succeeds on 7 of them. 

8 Conclusion 

In this paper, we present an incremental algorithm for computing CADs. A key part of 
the algorithm is an Intersect operation for refining a given complex cylindrical tree. If this 
operation is supplied with an equational constraint, it only computes a partial cylindrical 
tree, which provides an automatic solution for propagating equational constraints. We have 
implemented our algorithm in Maple. The experimentation shows that the new algorithm 
is much more efficient than our previous recursive algorithm. We also compared our im- 
plementation with the software packages Qepcad and Mathematica. For many examples, 
our implementation outperforms the other two. This incremental algorithm can support 
quantifier elimination. We will present this work in a future paper. 
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